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ATMOSPHERIC TEMPERATURE SOUNDINGS 
BY RADIOMETERS ABOARD SATELLITES 

by 

Roger Opstbaum 


Work has continued oh a survey of the methods for sounding the atmospheric 
temperature profile by remote measurements. The mainstream of this survey deals 
with infrared inversion techniques. However, the emphasis for this six-month 
period has been placed in reviewing related topics: sounding in the microwave 

region of the spectrum, sounding in cloudy atmospheres, and measuring sea-surface 
temperatures remotely. This survey was suggested by Dr. Joseph S. Hogan to pro- 
vide scientists at the Goddard Institute for Space Studies, who are developing a 
remote- sensing program for use with their dynamical models, with a compendium of 
previous research papers in the remote sensing field. 

At the end of the previous six-month period, the survey, which had proceeded 
chronologically, had gotten as far as Wark and Fleming's (1966) empirical orthogo- 
nal method (explained in the previous report of this series). This appeared to be 
the most recent non-iterative approach to the problem. A number of papers in 
recent years have chosen a method along the lines of Chahine (1968). The Planck 
function B[v,T(p)] found in the radiative transfer equation 

I(v,p) = e B [v,T(p q )] T(v,p o ) + J B [v,T(p) dp (1) 


in which I is radiation received at the satellite for wavenumber interval v at 
pressure p, e is the emissivity of the lower boundary which is at pressure p Q , 

T is the absolute temperature and T is the transmissivity, was often linearized 
in previous inversion schemes. Chahine, however, made use of the fact that the 
dependence of the transmission on the temperature was negligible in comparison to 
the dependence of the Planck function on temperature, and rather than seeking a 
linear form for equation (1), he was able to put it into the form 


I(v,p) 


X [T(p o )] A(v,p o ) + 


PP 

X[T(p)]K( v ,p) dp 

Jp o 


( 2 ) 


with "factors" X, A, and K. The above equation was solved for X(T) in simulation 
studies by creating a temperature profile T(p), determining its associated radi- 
ance profile, and than using the following iterative procedure to retrieve the 
temperature profile (quotes directly from Chahine, 1968): 



2 


1. "Select a set of ten sounding frequencies corresponding to ten pressure 
levels." Hogan (1972) advises me that numerical experiments have shown that due 
to the overlapping of the atmospheric weighting functions for carbon dioxide, in 
practice only 5 or 6 independent sounding frequencies can be chosen. 

2. "Make an initial guess T # (p) and solve equation 2; check the residuals 
and (R n ) av, and satisfy convergence." 

3. "Compare the results with the exact answer T(p)." The results are 
interesting. The guess converges quite rapidly to a rather accurate solution. 
Also, each isothermal guess converges to the same solution after about 5 itera- 
tions. However, the uniqueness of this solution, or, for that matter, any solu- 
tion, is still questionable, considering the nature of the weighting functions. 

The iterations have the effect of producing solutions which are accurate to only 
a few degrees from the actual profile. 

Several papers dealing with microwave sensing have been included in the sur- 
vey. Among them was the paper considered to be the pioneer paper in the field 
(Meeks and Lilley, 1963). This paper discussed results of calculations of opacity 
due to atmospheric oxygen in millimeter wavelengths and emphasized the complica- 
tion of the radiative transfer in this region of the spectrum due to Zeeman split- 
ting caused by the influence of the Earth^s magnetic field. The paper also out- 
lined a program to make the temperature measurements in the microwave region. 
Meanwhile, S tael in (1969) reviewed the work and applications of remote sensing in 
the microwave region. Among other important points made, he stated "microwave 
sensors provide the only remote sensing technique capable of measuring atmospheric 
temperature profiles in the presence of clouds." In view of the difficulty clouds 
present to people working in the field of infrared remote sensing, microwave tem- 
perature sensing will have to be given more emphasis and research than it has been 
given in the past. 

This is clearly seen to be true when a paper is encountered which discusses 
clouds and their effects on infrared sensing. For example. House (1968) has tried 
to introduce clouds into King's (1964) nonlinear method* He concludes— that accu- 
rate temperature soundings can only be made down to cloud-^top level,\a fact well 
established. before his study. \ 

Recognizing that clouds contaminate infrared data. Smith et al. (1970) have 
proposed a method for deleting cloudy data from a mixed data set for the 'purpose 
of making sea surface temperature measurements. The method can be termed 'a. sta- 
tistical-histogram method and it makes use of the fact that the brightness tem- 
perature for the clear atmosphere modal peak radiance is, statistically, the most 
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probable surface temperature# By collecting satellite data for a small latitudi- 
nal-longitudinal patch over a short period, a statistically random sample can be 
collected, which may or may not contain perfectly clear or cloudy data and prob- 
ably contains some mixture of each, so that a "clear mode" can be established. The 
method was successful enough to be considered for operational use. 

The survey is expected to be completed in the early fall^ 


ORBITAL DETECTION OF STRATOSPHERIC AEROSOLS 

by 

Gregory L. Matloff 

Research currently underway includes a feasibility study of the orbital 
detection of aerosols in the lower stratosphere. This research is under the 
direction of Dr. J. Hansen of the Goddard Institute for Space Studies. It 
appears that such effects as Rayleigh (molecular) scattering and refraction can 
be accounted for and that occultation of Stars by the Earth's atmosphere can be 
observed from satellites to determine local aerosol concentrations. 

MONTHLY PRECIPITATION CHARTS FOR THE WORLD 

by 

Roger Opstbaum 

Mean monthly precipitation charts have now been completed for all 12 months 
of the year. The data have been amassed from 6 sources (U. S. Department of Com- 
merce, 196 9a j U.S. Department of Commerce, 1969b; Critchfield, 1966; Miller and 
Thompson, 1970; Arakawa, 1969; and Wallen, 1970) and are plotted on a grid display 
according to latitude and longitude. 

As has been discussed in previous reports, because of the lack of precipita- 
tion data for the ocean areas, these charts cannot be used for such purposes as a 
heat budget study for the Earth. However, the great density of data over land 
areas makes the precipitation study purposeful and desirable. Its potential 
utility has been discussed in previous reports. 

The results confirm our previous general notions about global precipitation 
patterns. Furthermore, their detail over land areas and method of display shed 
some light on the mechanisms involved in seasonal and orographic causes of pre- 
cipitation. 



4 


Global precipitation patterns can probably best be described as patchy, with 
the greatest areas of high precipitation located in the subtropics and tropics* 
These are semi- permanent regions and they follow the movement of the Inter-Tropical 
Convergence Zone rather closely, although not all of them are located at the ITCZ. 
These patterns are responsible for Tropical Rain Forest type climates. 

Other precipitation types which are clearly illustrated in these charts are 
monsoonal (seasonal) and orographic precipitation mechanisms. The most striking 
example of a combination of both types of mechanisms can be found in the region 
around Cherrapunji, India. During the July monsoon, an average of greater than 
2000 mm is recorded for that area, while in some winter months, the amount of pre- 
cipitation recorded is less than that recorded in the eastern United States. 

More permanent orographically-induced precipitation is found, for example, 
in the Pacific Northwest in the region of the Cascades where precipitation values 
persistently higher than the surrounding areas are recorded the year found. This 
annual condition is probably due to the prevailing westerly winds from the ocean 
in that region. 

Most of the higher latitudes are not influenced by special topographic fea- 
tures and have a surprisingly constant amount of precipitation throughout the 
year, a condition more prevalent in the Northern Hemisphere than in the Southern 
Hemisphere. In the Southern Hemisphere, the transition from tropical to polar 
conditions is much more pronounced, possibly due to the large oceanic area, and 
this is reflected in precipitation amounts by a transition from heavy equatorial 
precipitation to the rather dry Antarctic continent. 

There are also numerous semi-permanent and permanent desert areas. Many sta- 
tions in the permanent deserts have an average of 0 mm for the entire year. It 
is rare, even in the cases of the semi-permanent deserts, for more than 50 mm 
precipitation to be recorded for one particular month in any desert region. 

DETERMINATION OF THE CLOUD STRUCTURE OF A PLANET 
FRCM REMOTE POLARIZATION MEASUREMENTS 
by 

Roger Opstbaum 

A study has begun in which the polarization of radiation arriving at Earth 
from selected points on a distant planet is calculated. These calculations are 
performed for various phase angles, optical depths, and models of the distribution 
of the multiple scattering particles. These multiple scattering computations 
involve the doubling method (van de Hulst, 1963; Hansen, 1971) together with the 
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method of successive orders of scattering (Uesugi and Irvine, 1970) which are 
being used by Dr. James E. Hansen at the Goddard Institute for Space Studies, who 
is collaborating on this study* 

Light undergoing multiple scattering by spherical particles becomes polarized 
in such a way that the polarization is a direct indication of the microstructure 
of the particles. Hence, if we want to determine the size distribution of the 
reflecting particles in a planetary atmosphere from the polarization characteris- 
tics of reflected sunlight, we can assume that the reflecting surface of the atmo- 
sphere under consideration is composed of spherical particles. It has been shown 
that the Venus clouds are composed of spherical particles (Hansen and Arking, 1971) 
This study can be applied most directly to Venus, where the cloud properties are 
still unresolved. Thus, by trying out several different simple models for types 
of scattering (e.g., Rayleigh, Mie, or some combination of the two for different 
optical depths) taking place in the Venusian clouds, and by comparing the result- 
ant theoretical polarization with observations of polarization for Venus (Coffeen 
and Gehrels, 1969), some insight can be gained as to the physical nature of the 
Venusian clouds. Also, by studying the latitudinal- longitudinal variation of the 
polarization, a better picture of the horizontal inhomogenities of the Venus 
clouds may become available. 

So far, intensity and polarization have been calculated for the cases of 
pure Rayleigh and pure Mie scattering. Currently, contour plots of these quanti- 
ties are being constructed. The intensity values for Rayleigh scattering have 
been found to be in total agreement with an independent source (Sekera and Viezee, 
1961). 


A MICRGDENSITCMETER ANALYSIS OF MARINER 6 AND 7 MULTICOLOR 
PHOT (METRIC PHOTOGRAPHS OF MARS 
by 

Gregory L. Matloff 

(Abstract of a paper submitted to Astronomical journal .) 

A photometric analysis of Mariner 6 and 7 three-color (B, G, R) photographs 
of Meridiani Sinus, the South Polar region, and Hellas, is performed using a 
microdensitometer. Hellas is apparently more specular in nature than other 
Martian terrain features and the dark areas near the Martian South Polar cap 
appear much darker than dark regions near the equator. Although much of the 
difference is due to AGC effects, some of the difference may be due to the 
different natures of the two terrain types. Two hundred and thirty-three points 
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in Meridiani Sinus are mapped in at least one color and nine points are mapped 
•in three colors. There is an apparent direct correlation between R/G reflect- 
ance ratio and height in the Meridiani Sinus region, based upon a comparison 
of R/G for four points ih Meridiani Sinus and earth-based radar measurements of 
that region. 

TECHNIQUE FOR PHOTOMETRIC DETECTION OF EXTRASOLAR PLANETS 

by 

•* 

Gregory L. Matloff and Alphonsus J. Fennelly 

(Abstract of a paper presented at a meeting of the American 
Optical Society in April 1972.) 

O'Neill (1968) has discussed the possibility of building an interferometer 
array of mirrors in space having an optical resolution equivalent to that of a 
5000 in. telescope. The applicability of such apparatus to detection and class- 
ification of extrasolar planets, using UBV photometry, has been discussed 
(Matloff, 1971). Utilizing Roman's (1959) estimate of the apparent magnitude of 
hypothetical planets circling Alpha-Centauri and O'Neill's discussion of the 
primary star's diffraction pattern, the S/N of the planet observed in the star's 
diffraction pattern and the photons received from the planet per unit time can 
be estimated. For an effective telescope diameter of 5000 in., S/N is about 0.1 
for earthlike planets circling Alpha-Centauri and 0.01 for Tau-Ceti or Epsilon- 
Eridani. With ah array of three widely separated 150 in* mirrors, about 100 
photons sec ^ are received in the case of the planet circling Alpha-Centauri A 
in the V and B photometric channels and about 70 in the U. In the case of Tau- 
Ceti, or Epsilon-Eridani, the planet flux is about an order of magnitude lower, 
and 70 percent lower for Alpha-Centauri B. Although O'Neill's discussion was 
concerned with mirror arrays in space, it is possible to construct an earth-based 
optical interferometer in an arrangement similar to that used in radio astronomy 
(Swenson, 1969), using optical instead of electronic methods to correct for sig- 
nal phase lag. If techniques of monitoring variations in the seeing can be 
developed, it may be possible to build an array of this type on a high terrestrial 
mountain instead of in space. 

American Elsevier Publishing Company, New York, New York. 
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ASYNCHRONOUS METEOROLOGICAL DATA INSERTION BY MEANS 
OF A COORDINATE TRANSFORMATION 

by 

Isidore M. Halberstam 

Much progress was made during the past six months in testing the method of 
data initialization by using continuous asynchronous data only. To make this 
report comprehensible as an individual unit; a brief review of the problem and 
the proposed method of solution are necessary. 

Introduction 

In I960, the Air Weather Service published a manual describing how weather 
analysts could best deal with reports submitted by air reconnaissance at times 
which were not simultaneous with standard synoptic reports^ The advised pro- 
cedure was to advect the data point, as well as possible, up to the time of 
synoptic reporting and to redraw the synoptic chart with the aid of the new data 
point. If a change in one of the readings (e.g., pressure) was indicated by 
changes in observed neighboring systems, the alteration was to be taken into 
account during advection. This procedure, as one could easily surmise, is 
highly subjective and is useful only if a small number of data points are avail- 
able. But since large bodies of asynchronous data were not readily submitted 
at that time, these shortcomings were not considered crucial. 

As the decade progressed, the problem of asynchronous data insertion took 
on a new meaning. With the advent of both satellite observations and numerical 
models of the atmosphere, the problem became one which involved large amounts of 
precious information which was needed as initial conditions for the numerical 
calculations. Satellites could cover areas which previously were devoid of ob- 
servers such as the vast stretches of the South Pacific far from the shipping 
lanes. But they could not supply all their information at one time. Instead, 
they were able only to scan a comparatively small area at each point in their 
orbit. By the time they reached half-way around the earth, the data they had 
previously supplied was outdated and could not correctly be used in conjunction 
with later- supplied data. The numerical models, however, can only make use of 
data supplied for one particular time. It uses this information to produce a 
forecast. Any information supplied after the initial moment can no longer be 
used in the initial conditions, but will confront values produced by the model 
in its forecast. If the later information is similar to the forecast value, 
there is obviously no need for it, since the model apparently simulated nature 
well enough. If the observation is vastly different from the forecast, the model 
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has obviously done poorly, due perhaps to a poor knowledge of the initial field. 

One feels intuitively that the later- supplied information could help the model's 
forecast. Just how to make use of it most efficiently has produced recent research 
activity. 

The most obvious procedure to handle asynchronous data is probably the 
simplest. It calls for simply inserting the data as it is received into its 
proper location. At each time-setp of the model, the satellite supplies data at 
a different grid point. If this data is forced into that grid point, the fore- 
cast will be altered. It was feared, however, that such a technique could not 
work if the neighboring grid points had forecast values which were different from 
the injected data. In that case, the detractors thought, an imbalance, or shock, 
would result and destroy the forecast. Jastrow and Halem (1970) attempted to 
demonstrate that such fears were unfounded. Their imaginative experiment showed 
that, to the contrary, inserted data would bring the computed values closer to 
the txrue solution and would cause no undue instabilities. Their experiment con- 
sisted of allowing the Mintz-Arakawa UCLA two- layer model to produce a forecast 
for a large number of days. The model was started again with the same conditions 
but after 119 days the forecast was interrupted and the entire field sprinkled 
with errors, both random and systematic. This caused the ensuing forecast to 
depart drastically from the original forecast. This departure was measured by 
the RMS error over the entire field of velocities. The model was allowed to run 
once more with the perturbations introduced at 119 days, but values from the 
original forecast were introduced at each time step at appropriate grid points. 

The RMS differences between the original and perturbed fields rapidly decreased 
to acceptable levels, even when the initial perturbation was extremely large. This 
would seem to indicate that raw satellite data could safely be inserted without 
too much difficulty into an ongoing numerical model and could substantially improve 
the forecast. 

Morel, Lefevre, and Rabreau (1971), among others, found that Jastrow and 
Halem' s results may not be universally applicable. They attributed the favorable 
results to the use of only one model in the experiment. The values used to simu- 
late satellite data were values obtained from the same model during a previous run. 
Since the system W3s naturally attuned to these values, it could easily be guided 
back to producing values close to the original ones. If real data were to be 
inserted, or if a simpler model without too many degrees of freedom were employed, 
instabilities may very well occur. Besides, the adjustment period made short 
range forecasts turn out rather poorly, and a forecaster desiring a 12- or 24-hr 
forecast could net take advantage of the satellite data. To overcome these 
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disadvantages. Morel proposed an iterative scheme to perfect a 24- or 48-hr baro- 
tropic forecast. This involves inserting observed data at its appropriate time 
and space coordinates until the desired termination time is reached. At this point 
the program is reversed starting from the 24- or 48-hr forecast field and hindcast- 
ing to time 0 while again inserting the same data at its appropriate location. At 
time 0 the program is moved forward again) starting from the hindcast conditions, 
with data insertion taking place as before* The process should converge giving a 
desirable result after a few iterations* Despite its ef fectiveneSs> this method 
is severely limited from a practical standpoint* Aside from the philosophical 
problems involved in reversing supposedly irreversible physical processes^ it is 
quite an expense to iterate a numerical model of some sophistication even for a 
short term forecast. Morel’s test on a barotropifc model is about as sophisticated 
as one can get without running up high expenses for computer time. 

A series of other proposals for the assimilation of asynchronous data were 
presented at the Symposium on Four-dimensional Data Assimilation held at Princeton, 
N.J., in 1971. Some indicated that weighting factors could be used to ease the 
data into the model without creating unduly severe gradients or shock waves. Others 
contended that weighting the data does not reduce the HMS error to any substantial 
degree. A short survey of these presentations can be found in Kasahara (1972), and 
repetition here is unnecessary. 

Coordinate transformation 

All the above recommendations have one object in common. That is, they all 
suggest means to insert data into an already specified model. It is hoped that 
despite the poor initialization in the model a better forecast can be produced if 
accurate observations are fed consistently into the model. But this expectation 
may not be borne out in a large number of cases. When simpler models are used, 
the introduction of new data may create the feared shock waves by constructing 
sharp gradients in the vicinity of the insertion. In sophisticated models, smooth- 
ing mechanisms or pseudo-viscosity terms can help keep the unstable modes in check 
while allowing the influence of the inserted data to spread. The ability to cope 
with the instabilities becomes, then, highly dependent on the type of numerical 
method employed. 

Current research seems to be moving away from smoothing and pseudo-viscosity 
toward implicit methods of solution. What effect insertion may have on these 
schemes is unknown, and no general conclusions can be drawn no matter how many 
schemes are tested. The reason for this lies in the mathematics of insertion. 

From the point of view of mathematics, once sufficient initial and boundary condi- 
tions are specified, the solutions to the differential equations of the model are 
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determined. Any "new" information is extraneous and can only be accepted if that 
information is compatible with the determined solution* If the information is not 
commensurate with the solution, the field is over-specified and either some of the 
specifications must be rejected or the solution will not exist. The same holds true 
when inserting data into a numerical modeli Since initial conditions are given, 
later insertion can only be considered overspecification. The only reason investi- 
gators have found success in data insertion experiments is that they are working 
with numerical analogs rather fhan analytical equations. These allow the inserted 
data to usurp the role of the initial conditions and slowly replace the solution 
based on the original conditions with the solution derived from the inserted data. 
Thus, it is hoped that the model will reject the first specification and keep the 
specification due to the inserted data. Analytically^ however^ this is still over- 
specification and can lead to destruction of the solution. It is like trying to 
solve six equations for five unknowns. If no two equations are linearly dependent, 
no solution can exist. The only way to insure a solution is to remove one of the 
specifications. 

To overcome the problem of "overspecif ication", it may be worthwhile to initial 
ize the field in a manner different from traditional methods. All investigators 
admit that the initial field based on poor observations ought to be modified by 
later observations based on satellite measurements. But calculations cannot begin 
unless there is some initial field. So, although sparse data is available over 
large portions of the globe for any given time, interpolations and extrapolations 
are performed based on data in neighboring areas. The later, more reliable data 
must then overcome the initial poor specification before the forecast can be made 
more reliable. If, however, the initial data could be ignored, and the later data, 
asynchronous though it may be, serve as the new "initial" field, the forecast may 
then be more reliable to begin with. This may be accomplished by means of a co- 
ordinate transformation in time. The satellite data is asynchronous, i.e. , it is 
a function of space and time. It can be made synchronous if the time coordinate, t 
is exchanged for a new coordinate, T , which is also dependent on space and time. 

It is constructed to coincide with the satellite path, so that all the data is fur- 
nished at the line t = 0 and is hence "initial". This coordinate is completely 
analogous to the <r-coordinate invented by atmospheric model builders to compensate 
for topographic features which could affect the traditional vertical coordinates 
z and p. The T-coordinate should accomplish the same task by smoothing out the 
"bumps" in the initial data curve. Thus, if data is given along a curve 
G(t,x,y) = 0, or, alternatively t =■ T(x,y), we may devise a time coordinate 
t = t - T(x,y) (or t/T(x,y) so that t o o (or 1) along the data curve, and the 
curve can indeed serve as initial data in the new coordinate system. 



11 


The coordinate transformation will, of course, alter the equations. If 
numerical means of solution are employed, it may be necessary to alter the scheme 
to fit the new terms. Truncation errors may increase if additional derivatives 
are necessary. On the other hand, stability requirements may be relaxed if coeffi- 
cients appear that can have the effect of shortening the time step. 

As an example, consider the equation: 


9P 

dx 


8£ + P 

9t 


in dimensionless variables, with initial condition P(x,0) “ 4x [exp (,5x)]. The 
solution here is P(x>t) « 4(x + t) exp [.5 (x - t)]* If information were not 
available at P(x,0) but along some curve P(x, t(x)), where t(x) *» 2>/x, say, we 
could still determine a solution to the equation. If we transform coordinates, 
such that (x, t) is replaced by (£,t), where £ " x and t » t - 2v/x, we get a 
new equation in terms of £ and r: 

3p 0P f 1 + -n/1 ^ „ 


Now, however, the condition that was specified along P(x, t(x)) is a condition 
at P(£,0) which appears as an ordinary initial condition in (£, T ) coordinates. 

The solution obtained should still be the same but will be in terms of (£, T ) 
instead of (x, t) . Of course, the transformation can be used to express the solu- 
tion in terms of (x,t) if necessary. 

An experiment was performed using the numerical analog of this equation. 
Results of the experiment were described in the previous report submitted in 
December 1971. 

The Wave equation 

A similar experiment was conducted with the wave equation and with a one- 
dimensional gravity wave equation. The wave equation was first used in its one- 
dimensional form and then as a two-dimensional equation. The results proved simi- 
lar in both cases, so only the two-dimensional equation will be discussed here. 

The wave equation is in the well-known form 


2 2 2 
2 a (b 2 9 <t> 9 <t> 

c 2^* c 2° 2 * 
9x 9y 9t 


(la) 


where c represents the phase speed of the wave. Its form is well-known, and the 
scheme is conditionally stable. If i refers to the time index, J, the x-index, 
and f, the y-index, the numerical equation is: 
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The exact solution to the equation depends on the Initial and/or boundary condi- 
tions. In this experiment, the boundary and initial conditions were chosen to 
coincide with a solution equal to 

<J>(x,y,t) “ 3 sin .01 (x - ct) + 3 sin .01 (y - ct) (2) 


The transformed equation looks rather complicated because of the second 
order derivatives. The transformation was given as § «=> x, tj “ y, and 
T = t - 2S Ik + y + 15. In the (£,t), T ) coordinate system the equation becomes: 



Some difficulty was encountered trying to find the numerical analog to this 
equation. Unlike the first-order equation, this equation has cross derivatives 
between time and space. To use purely centered differences would somehow involve 
using an implicit scheme since there would be more than one variable at the i + 1 
time step. This is technically feasible but the increased computation time and 
the difficulty of the programming would make the original equation much more 
appealing than the transformation. To avoid this, a backward difference was used 
on the cross-terms so that there was only one variable to be solved at i + 1. 

But this also increased the truncation error considerably. However, when a 
second iteration was performed using the values generated by the backward differ- 
encing and employing a strictly centered differencing, the truncation error was 
reduced considerably. The numerical equation for the first iteration was then: 
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1+1 


. 1+1 


and 


The m\ly unknown quantity here is B » All others, including $ n 
1+1 Jj* 3 x > Jt 

, are known from previous time steps. The second iteration affects only 
Jj*~ 1 i+1 i-1 

the third and fourth terms on the left hand side where 0. . and 4> . receive 

the indeces j+1 and £+1. 

To compare results, the transformation was used to transform values from 
different T-times to t-times. Since values were available only at finite steps 
in T most values obtained were only close to, but not exactly equal to, the 
real time t. However, since the time step used was rather small, this should 
not be noted as a significant cause for departure. Parameters used in the cal- 
culation were c = .8, At * .9, Ax « Ay = 2.1. The third part of the experi- 
ment involved data insertion after perturbing the initial field with random 
errors equal to £ 1/8 the values of the initial field. The insertion was done 
similarly to the manner discussed in the previous report. The results in terms 
of average absolute truncation errors for a few time steps are shown in Table I. 
(Truncation error here is simply the difference between the exact solution and 
the numerical solution). 

As can be seen, the centered difference method using exact initial condi- 
tions allows for the lowest truncation errors with a very slow growth rate of 
these errors. The results using equation (3b) gives truncation errors two orders 
of magnitude larger than the results using equation (lb). This is due to three 
factors: First, as mentioned before, the backward differencing causes an increase 

in truncation error; second, there are more terms with derivatives adding to the 
total; third, the coefficient of the time derivative on the right-hand side of 
(3b) magnifies the time step AT. This term is less than one, which is the same 
as multiplying At by a number more than one, which will increase the truncation 



14 


Time 

Table Ii Average Absolute 

Errors for <}> 

Data insertion 

Centered 

differencing 

< lb > xlO- 6 

A sy nchronou s , 
no iteration 

< 3b > xlO' 4 

Asynchronous, 
one iteration 

x 10" 6 

49.50 

3.0363 

3.6101 

5.8603 

20.207 

50.40 

3.0820 

3.6966 

6.0537 

20.635 

51.30 

3.1255 

3.7824 

6.2489 

21.041 

52.20 

3.1668 

3.8673 

6.4454 

21.422 

53.10 

3.2057 

3.9513 

6.6435 

21.780 

54.00 

3.2423 

4.0344 

6.8426 

22.117 

54.90 

3.2764 

5.1164 

7.0425 

22.432 

55.80 

3.3080 

4.1974 

7.2432 

22.724 

56.70 

3.3371 

4.2773 

7.4444 

22.986 

57.60 

3.3636 

4.3559 

7.6461 

23.222 

58.50 

3.3874 

4.4333 

7.8481 

23.439 

59.40 

3.4085 

4.5094 

8.0502 

23.633 

60.30 

3.4269 

4.5841 

8.2524 

23.794 

61.20 

3.4425 

4.6574 

8.4544 

23.934 

62.10 

3.4553 

4.7292 

8.6561 

24.056 

63.00 

3.4652 

4.7995 

0.8573 

24.159 

63.90 

3.4722 

4.8682 

9.0578 

24.226 

64.80 

3.4763 

4.9352 

9.2575 

24.274 

65.70 

3.4774 

5.0006 

9.4564 

24.317 

66.60 

3.4756 

5.0642 

9.6543 

24.335 


error while requiring stricter linear instability criteria. The iteration proce- 
dure reduces the errors to the same order of magnitude as equation (lb), but the 
rate of error growth is still larger. For the data insertion method, the shock 
effect is quite evident since the errors are from three to four times the size of 
the variable. The solutions obtained without inserting data, in fact, were better 
behaved and their departure was not as great. This seems to indicate that at 
least in simple models, without any sophisticated damping devices, the threat of 
shock waves being created is quite real when raw data is inserted. Using the 
transformed coordinates, however, circumvents this problem. 
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One-dimensional free surface gravity waves 

A test experiment was performed on primitive barotropic equations whose 
variables depended only on x. The linear counterpart of this model is dis- 
cussed by Kurihara (1965). The equations used here are: 


3u 

3t 


+ mu 


9u 

3x 


3<j> 

£v • ^ 


(4a) 


3v 

at 


+ 


mu 


8v 

3x 


fu 


(4b) 


—• + mu ® vfU - m(gH + <t>) §~ ( Ac ) 

where u is the westerly component of velocity, v is the southerly component, 
f is the Coriolis parameter, U is the average westerly component (16 m sec ), 
H(y) is the mean height of the free surface (5800 - “ y), m is the map factor, 
and $ is the deviation of geopotential from the mean, m and f depended only 
on the y-ditection while u, 4> and v only varied in the x-direction. The numeri- 
cal scheme employed was the same employed by Shuman (see Haltiner, 1971, p. 227) 
with simplifications made for the one-dimensional case. We shall use Shuman's 
notation which, briefly summarized, gives: 

%<A J4% + A * S £ (A j4% ‘ A J-V 


and so 


A ™ “ ^< A j+ i + + A j-i>^ A x - ite ( Vi 

The numerical analog of (4) is given by 

> ■ — — 5C 

-t 


A J-1* 


-X-X . TX-X 
u «»-muu + r v - 4> 
t K X 


~t -X-X 7X-X 

v. « - m u v - iu 

t X 


(5a) 

(5b) 


o - m X u X <t> + v 54 ? 1 ^ - tn X (gH + £ X )u (5c) 

t X X 

The second term on the right of (5c) is equivalent to the first term on 

3H 

the right in (4c) and is derived from the equality vg ^ <= - fvU which defines 
the dependency of H on y. The averaging symbol over m and f is academic in 
this model, since they do not vary in the x-direction. Using a time step of five 
minutes and Ax * Ay ° 200 km, forecasts were made for 12, 24, 36, and 48 hr, 
and the heights H + — were plotted by the S-D 4060 plotter at the Goddard 

s 
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Institute for Space Studies (GISS) in New York City. Extra smoothing was inserted 
every 50 time steps to keep the curves smooth and to avoid any instabilities that 
might arise. The initial conditions were simple sinusoidal waves for all three 
variables with their respective tanges as 12 ^ u ^ 20, — jj“ x 9.8 ^ v ^ — g x 9.8, 

- 980.0 <t> 980.0. The boundaries in the X-direction were cyclic so that the 

last grid points always coincided with the first. 

The experiment was conducted in the same manner as the previous ones, but 
since no exact solution was available here, the results of equation (5) were used 
as the true solution and as data for the initialization of the asynchronous equa- 
tions and for insertion during the third part of the experiment. 

For the asynchronous experiment, the transformation was given in analytical 
terms. The transformation was 

t <= t + 1 - 10 3 (x + y) 

This coincides with a satellite speed of approximately 1 km sec 1 . The transformed 
equations were then: 



[*<-mu ff + fv - §f ) - T x (fUV “ mU tf " m(gH + 40 if } 3/0 


(7a) 


3<j) 

w 



fu]/T 


(7b) 


3<t> P • d (b . 9u. 

« [ T (-mu g|- + fUv - m(gH + 0) fig) 

- ^(f (gH + 0) mv - m 2 u(gH + 0) ^ - m(gH + 0) ^3/D (7c) 

.22 • 8t 9t 

where D = T - mT^gH +0), T = 1 + mu 3^> T x =: 8^ 

The numerical counterpart of (7) is exactly similar to (5), since T and 
D were calculated separately and are not involved in the smoothing of the products. 
Computation was carried out until all grid points had been calculated up to 48 hr 
of real time. An interpolation was used to recover the fields of u, v, and 0 
at 12, 24, 36, and 48 hr of real time. The interpolation consists of using equa- 
tion (6) to test at each T-step if the value of T coincided with t = 12 hr, say. 
When that value of t was located between two values of T , a linear interpola" 
tion was carried out which estimated the asynchronous value of the variable for 
the required time. If larger time steps are used, or if the fields undergo quick 
change, better methods of interpolation may be necessary. The same method was 
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employed (in reverse) to derive the field at T * 0, when the field at all t 
was known. 

The boundaries in the asynchronous case gave another sort of problem. Since 
the boundaries were assumed cyclic for the synchronous case, it was desired to 
keep them cyclic in the asynchronous case. However, for any given T-step it was 
no longer true that the last grid points have the same value as the first ones, 
since, by equation (6) they were not at the same real time. It was easy enough to 
define the "later" boundaries from earlier ones, but the reverse was virtually 
impossible unless some sort of implicit scheme or iterative method was performed. 

In practice^ the boundaries were determined by interpolating from the synchronous 
results at each time step. This did not upset the solution too badly, nor did it 
nudge the results closer to the true ones to any degree. 

The third part of the experiment involved perturbing the initial field by 
random errors which, in magnitude, were 25 percent of the velocities and ^ 100 m 
for the height field. Data was inserted at each time step according to equation (6) 
by setting T «= 0. The set of grid points at each time step that coincided with 
the x and y values derived were fed interpolated data. This meant that at each 
time step, a diagonal swath of grid points received information until the entire 
field was covered in about 2,5 hr. The results of the experiment are shown in 
Figures 1-13. The choppiness in the asynchronous graphs is due mainly to the inter- 
polation. Otherwise, the angles and troughs are accurately depicted at all times. 
The data insertion, as expected, does poorly in the early forecasts but improves 
toward the end of the period, although no new data is inserted after the initial 
2.5 hr. This is due in part to the smoothing which serves the dual purpose of 
cancelling out the random errors and spreading the influence of the inserted data. 

It is also due to the nature of the fixed slope in the y-direction which tends to 
act as a forcing function in equation (4c) so that by 48 hours the data insertion 
method gives nearly the same results as the error-free model. 

For a numerical comparison, values for u are given in Table II at grid points 
f a 7, j as 19-21, for times 0, 12, and 36 hr. Notice how once again the insertion 
method does better for the long-range forecast than for short-range. The effect 
of smoothing in cancelling the random errors can be seen by the improvement of the 
results where there was no insertion whatsoever. 

Note on smoothing 

Special mention must be made regarding the effect of smoothing on the differ- 
ent techniques. By smoothing properly, it is hoped, the high frequency computa- 
tional modes can be damped while the important long-wave physical mode is allowed 
to survive. Here, the same smoothing algorithm was used for both the synchronous 
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Table II. Values of 

u(m sec ^) 

for different 

times 

at 

L = 7, 

19 ^ J « 21 



Time (hr) 


0 

12 

36 

j = 19 

12.763 

9.320 

-10.442 

Synchronous 

20 

13.323 

9.320 

-10.442 

(no error) 

21 

14.000 

9.491 

-11.096 


19 

- - - 

9.191 

-10.295 

A synchronou s 

20 

- - - 

9.181 

-10.733 


21 

- - - 

9.895 

-11.023 

Initial 

19 

13.777 

7.585 

-10.056 

Random error 

20 

13.254 

7.534 

- 9.653 

(no insertion) 

21 

11.943 

8.396 

-10.080 


19 

13.777 

7.622 

-10.451 

Data insertion 

20 

13.254 

8.470 

-10.508 


21 

11.943 

9.312 

-11.111 


and asynchronous case. It was originally feared that such a procedure would 
create serious deviations because different data would be incorporated in each 
case. In the synchronous case, smoothing balances the value at a grid point with 
the values found in neighboring points at a given time. In the asynchronous case, 
the values in neighboring grid points are not values given at the same real time. 
Therefore, smoothing in the asynchronous case means smoothing in time as well as 
space. When excess smoothing was induced, the values began to diverge, in fact. 
The contours appeared flattened after 12 hr in both cases, but the values of the 
contours were significantly different. With the reduction of extra smoothing to 
only once every 50 time steps, the correlation improved drastically as can be 
seen from the plots. In the case of the insertion method on the other hand, in- 
creased smoothing brought the results in line much more quickly, for the reasons 
previously mentioned. Thus smoothing seems to be a damping technique which is 
beneficial to the insertion of data technique but detrimental to the coordinate 
transformation method. 
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Conclusions and proposals 

The results so far in this investigation have seemed very promising and call 
for future research. The work on the first order equation and the wave equation 
prove that shock waves are a threat in data insertion. The free surface equations 
show that data insertion is not useful for short-range forecasts, while the trans- 
formation procedure shows definite promise of being able to handle prediction 
from the time the entire field is covered. 

The next step in this research will be to undertake, as much as possible, a 
simulated procedural run using the asynchronous transformation. The two-dimen- 
sional barotropic equations will be used for a small region on the Shuman grid. 
Initial conditions, boundary conditions, and satellite data will be gleaned from 
the six-layer model of Shuman. The results will be analyzed both from a compari- 
son of 500-mb plots and from a comparison of error indicators such as the RMS 
errors 
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